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1,  Introduction 

This  report  is  the  final  report  of  the  Phillips  Laboratory  contract  F19628-93-K-0033.  This 
contract  was  funded  in  response  to  a  proposal  that  was  submitted  to  Phillips  Laboratory  in  July 
1993  with  the  title  Improved  Estimation  of  the  Earth’s  Gravity  Field.  The  starting  date  of  the 
contract  was  22  October  1993  with  an  ending  date  of  21  January  95.  On  May  31,  1994  a  request 
was  made  to  Anestis  Romaides  to  expand  the  scope  of  the  original  proposal  with  additional 
funding  requested.  In  July  1994  the  modification  of  the  contract  was  approved  and  the  contract 
ending  data  was  extended  to  22  April  1995. 

2.  Research  Activities 

The  proposal  submitted  in  July  1993  described  the  importance  of  improving  our  knowledge  of 
the  Earth’s  gravity  field  and  its  use  in  the  ultimate  development  of  a  world  height  system.  In  the 
context  of  the  original  proposal  the  Earth’s  gravity  field,  or  more  properly  the  Earth’s  gravitational 
potential  was  to  be  represented  by  a  set  of  potential  coefficients  and  techniques  for  the  improved 
estimation  of  these  coefficients  was  part  of  this  project’s  area  of  interest.  In  the  May  31,  1994  letter 
to  Mr.  Romaides  the  work  of  the  project  was  to  be  expanded  so  that  project  personnel  could 
participate  in  the  research  and  meetings  related  to  the  joint  gravity  and  geoid  improvement  project 
of  the  Defense  Mapping  Agency  and  the  NASA,  Goddard  Space  Flight  Center.  This  report  has 
been  written  to  describe  the  activities,  undertaken  by  project  personnel,  consistent  with  the  goals 
defined  by  the  proposals.  The  persons  involved  with  research  supported  by  this  project  are  Richard 
H.  Rapp,  Professor  Emeritus,  Christopher  Jekeli,  Associate  Professor,  Nagarajan 
Balasubramania,  Dru  Smith  Graduate  Research  Associates.  For  convenience,  this  final  report  has 
been  divided  into  sections  according  to  the  author. 

2.1  Research  Activities  of  N.  Balasubramania 

Mr.  Balasubramania  continued  his  activities  related  to  the  development  of  a  world  height 
system.  His  research  led  to  the  publication  of  Scientific  Report  No.  1,  in  April  1994,  "Definition 
and  Realization  of  a  Global  Vertical  Datum,"  PL-TR-94-2144,ADA283303,  the  abstract  of  this 
report  is  as  follows: 

Definition  and  Realization  of  A  Global  Vertical  Datum 
Nagarajan  Balasubramania 
Abstract 

The  requirements  for  defining  and  realizing  a  global  vertical  datum  (GVD)  with  a 
precision  on  the  dm-or  even  cm-level  has  become  essential  now,  due  to  the 
improved  measurement  accuracy  provided  by  the  space  geodetic  techniques.  This 
study  identifies  the  different  approaches  in  which  such  a  global  vertical  datum  could 
be  evolved. 

The  first  approach,  which  can  be  considered  as  an  ideal  approach,  uses  the  best 
available  geocentric  station  positions  from  present  and  future  space  geodetic 
networks,  highly  accurate  geopotential  model  and  surface  gravity  data  around  the 
space  stations  for  defining  and  realizing  a  global  vertical  datum.  The  other 
simplified  approach,  which  can  be  considered  more  of  an  operational  development 
with  a  short  time  framework,  is  mainly  dependent  on  the  widely  available  GPS  and 
DORIS  tracking  networks  and  accurate  geoid  height  models  for  the  GVD 
development. 
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After  a  review  of  the  various  heights  and  height  systems  that  are  in  use  today, 
detailed  mathematical  procedures  for  setting  up  the  observation  equations  to  define  a 
global  vertical  datum  in  a  least-squares  sense,  under  different  data  scenarios,  are 
.discussed.  A  test  computation  to  realize  the  first  iteration  global  venical  datum  with 
available  data  at  17  space  geodetic  stations  in  six  regional  vertical  datums  is  also 
included.  The  gravimetric  height  anomaly/undulation  computations  at  the  space 
geodetic  stations  required  in  the  test  computation  have  been  performed  using  both 
Modified  Stokes’  technique  and  least  squares  collocation  technique  combining 
surface  gravity  data  in  a  small  cap  around  the  stations  with  potential  coefficients 
from  OSU91A  model.  The  test  results  show  that  in  an  idealistic  approach  the  global 
vertical  datum  can  be  realized  to  an  accuracy  of  ±5cm  an  its  connection  to  the 
regional  vertical  datums  to  an  accuracy  of  ±5cm  to  ±23cm.  The  estimated  values  of 
separation  between  different  regional  vertical  datums  agree  mostly  well  with  the 
results  reported  by  various  geodesists  and  oceanographers  based  on  their  regional 
studies. 

The  Balasubramania  report  also  appeared  as  Report  427  of  the  Department  of  Geodetic  Science 
and  Surveying,  April  1994.  The  report  was  also  submitted  to  the  Graduate  School  of  The  Ohio 
State  University  in  partial  fulfillment  of  the  requirements  for.  the  Degree  Doctor  of  Philosophy. 
Basasubramania  received  his  PhD  from  Ohio  State  in  March  1994. 

2.2  Research  Activities  of  Richard  H.  Rapp 

During  this  time  period  Professor  Rapp  served  as  the  adviser  to  Balasubramania  and  carried  out 
activities  related  to  geoid  and  gravity  model  improvement. 

The  interest  in  vertical  datums  was  continued  in  this  project.  In  December  1993  made  two 
presentatioris  in  Wormley,  England,  related  to  project  interests.  The  first  presentation  was:  A 
Global  Vertical  Reference  Frame.  This  was  at  the  special  meeting  of  the  group  involved  with  the 
Geodetic  Fixing  of  Tide  Gauge  Bench  Marks.”  The  meeting  of  this  group  was  followed  by  the 
meeting  of  the  International  Association  of  Geodesy  Special  Study  Group  5.149,  Vertical  Datum 
Investigation.  At  this  meeting  Rapp  made  the  presentation:  A  World  Vertical  Datum  Proposal. 
Copies  of  the  viewgraphs  used  in  this  talk  are  given  in  the  Appendix  of  this  report. 

Upon  invitation,  Rapp  prepared  the  following  paper  for  presentation  at  the  meeting  for  the 
International  Symposium  on  Marine  Positioning,  Hannover,  Germany:  Vertical  Datums  in  Land 
and  Marine  Positioning.  This  symposium  was  held  in  September  1994.  The  complete  text  of  this 
paper  is  contained  in  the  Appendix  of  this  report.  The  paper  was  also  published  in  the  Proceedings 
of  the  Symposium  published  by  the  Marine  Technology  Society,  Washington,  D.C. 

Rapp  has  also  been  involved  in  the  joint  DMA/GSFC  project  for  gravity  field  and  geoid 
improvement.  He  is  the  chair  of  the  steering  committee  and  participates  in  the  activities  of  the 
working  groups  charged  with  the  implementation  of  the  procedures  that  will  lead  to  a  significantly 
improved  degree  360  model.  During  this  contract  period,  Rapp  chaired  two  meetings  (April  5  and 
November  1,  1994)  of  the  steering  group.  He  prepared  minutes  of  these  meetings  which  were 
distributed  to  all  parties  of  the  joint  effort.  In  addition  to  the  steering  group  meetings  Rapp  attended 
working  group  meetings  held  8-9  June  94  (NASA/GSFC)  and  1-2  August  94  (DMAAC).  At  the 
July  meeting  Rapp  made  a  presentation  on:  Terrain  Base  (NGDC)  and  Ice  Thickness  Information. 
This  informal  presentation  described  comparisons  made  with  several  different  elevation  files. 
Specific  comparisons  were  carried  out  by  Dm  Smith,  Graduate  Research  Associate,  in  a  l°xl°  celi 
(NW  corner,  (1)=20°,  ?l=204°),  in  the  Hawaiian  area.  5’x5’  elevations  were  computed  with  a  3”x3” 
DTED  file  provided  to  Ohio  several  years  ago.  These  mean  elevations  were  then  compared  to  the 
elevations  in  three  data  bases:  ET0P05U,  T  Base  (a  version)  and  T  Base  (p  version).  The  Terrain 
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Base  data  sets  were  provided  to  us  by  NGDC  for  evaluation.  Statistics,  as  complied  by  Smith,  are 
given  in  Table  1 


Table  1 

Comparison  of  5’x5’  Mean  Elevations  From  Three  Data  Sets  to 
Ground  Truth  Values  Estimated  from  3”x3”  Values 


Model 


ET0P05U 

T  Base  (a) 

T  Base  (p) 

No.  of  Comp 

89 

87 

87 

Mean  Diff. 

530m 

475m 

155m 

Std.  Dev. 

±916m 

±7 18m 

±434m 

RMS  Diff 

1058m 

861m 

461m 

Min  Diff 

-1217m 

-730m 

-730m 

Max  Diff 

2548m 

2936m 

1499m 

No.  of  Diff 

'  37 

18 

5 

where 

Id  l>100m 

The  conclusion  from  this  test  was  that  Terrain  Base  ((3)  was,  in  this  area,  a  significant 
improvement  over  the  older  ET0P05U  and  the  initial  version  of  Terrain  Base. 

Rapp  and  R.S.  Nerem  of  NASA/GSFC,  prepared  a  paper  for  the  Joint  Symposium  of  the 
International  Gravity  Commission  and  the  International  Geoid  Commission  that  was  held  in  Graz, 
Austria,  September  1 1-17,  1994.  The  title  and  abstract  of  this  paper  are  as  follows: 


A  JOINT  GSFC/DMA  PROJECT  FOR  IMPROVING  THE  MODEL  OF 
THE  EARTH'S  GRAVITATIONAL  FIELD 

Richard  H.  Rapp 

Department  of  Geodetic  Science  and  Surveying 
The  Ohio  State  University 
Columbus,  Ohio  43210 


R.  Steven  Nerem 

NASA  Goddard  Space  Flight  Center 
Space  Geodesy  Branch,  Code  926 
Greenbelt,  MD  20771 


ABSTRACT 

The  U.S.  Defense  Mapping  Agency  and  the  NASA  Goddard  Space  Flight  Center 
with  the  aid  of  other  organizations  such  as  The  Ohio  State  University  are 
cooperating  in  a  joint  effort  to  determine  a  significantly  improved  degree  360 
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spherical  harmonic  model  representing  the  Earth's  gravitational  potential.  This  new 
model  will  be  of  immediate  use  in  defining  a  geoid  undulation  model  that  will  be  the 
basis  for  an  enhanced  WGS84  geoid. 

The  development  of  the  new  model  is  driven,  in  part,  by  the  need  to  determine  an 
accurate  geoid  undulation  model  that  will  be  the  reference  surface  for  a  World 
Height  System  to  be  implemented  in  the  1996  time  period.  In  addition,  the  new 
geoid  estimation  will  help  satisfy  increasingly  important  studies  in  ocean  circulation 
(sea  surface  topography)  and  geodetic  positioning  through  GPS. 

The  new  model  estimation  will  incorporate  existing  and  new  satellite  data.  New 
data  will  include  GPS  tracking  of  Topex/Poseidon,  Doris  tracking  of  several 
satellites,  altimeter  data  from  Topex/Poseidon,  and  ERS-1  and  Doppler  data  from 
satellites  at  inclinations  not  covered  or  weakly  represented  in  previous  solutions. 

The  surface  gravity  data  to  be  used  in  the  new  solution  will  be  based  on  an 
updated  30'  mean  anomaly  data  base  developed  at  the  DMA  Aerospace  Center.  This 
new  data  set  will  incorporate  a  substantial  amount  of  new  data  that  has,  and  will, 
become  available  in  Europe,  the  FSU,  South  America,  Greenland,  Africa,  Asia  and 
Antarctica.  Anomaly  values  in  areas  such  as  Canada,  the  United  States  and 
Australia  will  be  based  on  updated  data  files.  This  data  will  be  used,  after  suitable 
corrections,  to  form  normal  equations  that  can  be  used  with  the  satellite  derived 
normal  equations. 

In  addition,  30'x30'  mean  anomalies  derived  from  the  Geosat  Geodetic  Mission 
satellite  altimeter  data  will  be  used  in  the  project.  The  file  will  be  merged  with  the 
files  based  on  the  surface  terrestrial  data.  In  areas  where  no  data  exists,  anomaly 
estimates  will  be  made  from  new  elevation  data  through  topographic  isostatic 
models  to  ultimately  yield  a  global  30'  anomaly  file.  In  addition,  an  updated  l°xl° 
anomaly  file  based  on  terrestrial  data  (both  land  and  ocean)  will  be  determined. 

The  final  stage  of  the  data  processing  will  be  the  development  of  several  degree 
360  models  using  different  data  sets  and  weighting  procedures.  The  current  plan  is 
to  use  existing  software,  for  the  combination  solution,  with  minimum  modifications 
to  assure  a  timely  effort.  Several  preliminary  models  will  be  made  available  to  the 
international  community  for  evaluation.  A  final  model  will  be  selected  based  on 
extensive  tests  of  the  preliminary  models.  The  final  model  and  accuracy  estimates 
will  be  released  in  mid  1996.  The  model  will  be  used  to  determine  accurate  geoid 
undulations  that  will  be  available  in  gridded  form. 

This  paper  will  be  published  in  the  proceedings  of  the  meeting  by  Springer  Verlag  in  the  lAG 
series.  Copies  of  the  complete  paper  were  provided  to  Phillips  Laboratory  and  DMA.  The  paper 
was  to  have  been  presented  by  Rapp  but  due  to  unforeseen  circumstances  he  was  not  able  to  attend 
the  meeting.  The  paper  was  delivered  by  Dr.  M.  Kumar  of  DMA. 

Additional  studies  were  earned  out  to  consider  precise  procedures  for  the  determination  of  geoid 
undulations  from  the  potential  coefficient  representation  of  the  Earth’s  gravitational  potential. 
Previous  discussion  on  this  evolved  knowing  the  difficulty  in  evaluating  the  disturbing  potential  at 
a  surface  (the  geoid)  interior  to  the  Earth.  An  alternative  suggestion  has  been  made  that  the 
disturbing  potential  be  evaluated  at  the  external  surface  of  the  Earth  and  then  downward  continued 
to  the  geoid  using  additional  information.  The  downward  continuation  can  be  done  using  equations 
developed  many  years  ago  by  Moritz  who  showed  the  relationship  between  the  geoid  undulation 
and  height  anomaly.  This  relationship  depends  on  Bouguer  gravity  anomalies  that  can  be  computed 
from  a  potential  coefficient  model  if  topographic  elevations  are  known.  Postulating  that  a  spherical 
harmonic  expansion  of  terrain  elevations  can  be  made  procedures  can  be  developed  to  represent  the 
correction  term  to  convert  height  anomalies  to  geoid  undulations.  In  a  simplified  case  correction 
terms  for  the  potential  coefficient  can  be  computed.  This  could  lead  to  the  determination  of  a 
modified  set  of  potential  coefficients  which  could  be  used  for  geoid  undulation  computation  using 
traditional  software.  The  problem  is  complicated  by  the  fact  that  evaluation  of  the  disturbing 
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potential  at  altitude  requires  knowledge  of  the  terrain  elevation.  This  significantly  reduces  the 
efficiency  of  computer  software  designed  for  the  fast  calculation  of  geoid  undulations,  from 
potential  coefficient  models,  on  a  grid.  The  best  way  to  carry  out  these  computations  has  not  yet 
been  determined. 

A  review  of  methods  for  the  calculation  of  geoid  undulations  from  potential  coefficient  models 
was  prepared  by  Rapp  for  this  contract  and  for  presentation  at  the  International  School  for  the 
Determination  and  Use  of  the  Geoid.  This  school  was  held  in  Milan,  Italy,  October  1994,  under 
the  auspices  of  Professor  Fernando  Sanso,  Director.  The  title  of  this  paper  is;  The  Use  of  Potential 
Coefficient  Models  in  Computing  Geoid  Undulations.  Copies  of  the  paper  have  been  provided  to 
Phillips  Laboratory  and  the  Defense  Mapping  Agency.  The  paper  will  be  published  in  the  bulletin 
series  of  the  International  Geoid  Service  in  1995. 

2.3  Research  Activities  of  Christopher  Jekeli 

The  traditional  method  of  spherical  harmonic  analysis  of  discrete,  global  ^avity  anomaly  data 
was  popularized  by  Colombo  who  reformulated  the  technique  for  FFT  utiliziation,  thus  offering 
the  ability  to  perform  very  high  degree  expansions  with  minimal  computational  effort.  However, 
the  method  is  still  based  on  a  discretization  of  the  integral  (a  quadrature  formula)  that  defines  the 
spherical  harmonic  coefficient  in  terms  of  the  anomaly  data.  Other  methods  based  on  least-squares 
collocation  and  least-squares  adjustment  (also  developed  by  Colombo)  similarly  discretized  the 
integral  (albeit  with  different  weights),  a  necessity  because  the  data  are  discrete. 

Each  of  these  methods,  then,  is  subject  to  an  error,  known  in  the  parlance  of  spectral  analysis, 
as  an  aliasing  error.  It  is  an  error  independent  of  the  noise  in  the  data,  but  does  depend  on  the 
amount  of  information  contained  in  the  data  at  frequencies  generally  higher  than  the  frequency 
defined  by  the  sampling  interval  (the  Nyquist  frequency).  In  theory,  and  if  certain  assumptions  of 
the  geo- statistical  nature  of  the  gravity  field  are  fulfilled,  the  least-squares  collocation  technique 
should  yield  the  best  (minimum  aliasing  error)  estimate  of  the  harmonic  coefficients.  The  aliasing 
error  can  be  estimated  empirically  using  simulation  studies;  and  some  simple  models  were 
proposed  for  the  OSU89  and  previous  models.  This,  in  turn,  led  to  simplified  weighting  schemes 
for  potentially  better  estimation  of  the  harmonic  coefficients. 

This  rather  ad  hoc  approach  to  harmonic  analysis  is  unsatisfactory  and  calls  for  a  more 
systematic  characterization  of  the  aliasing  error.  This  was  accomplished  by  Jekeli  who  derived 
explicit  formulas  for  the  aliasing  error  for  the  traditional  quadratures  formula,  and  by  analogy  for 
the  least-squares  adjustment  formula.  Also,  a  modified  estimation  was  proposed  that  reduces  the 
aliasing  error  by  determining  twice  the  number  of  harmonic  coefficients  (thus,  the  same  number  as 
number  of  observations)  at  only  slightly  greater  computational  cost. 

In  addition,  by  clearly  understanding  the  origin  of  the  aliasing  error  in  the  analysis  of  data  on  a 
sphere,  it  became  apparent  that  certain  types  of  data  smoothing  are  better  than  others  in  reducing 
the  aliasing  error.  Typically,  data  smoothing  is  the  result  of  uniformly  weighted  averaging  of 
gravity  anomalies  over  blocks  defined  by  latitude  and  longitude  lines.  The  aliasing  error  formulas 
show  how  this  procedure  fails  to  attenuate  a  part  of  the  high-degree  spectrum  that  then  aliases  the 
estimated  spectrum.  Averages  over  constant  spherical  caps,  though  spoiled  at  higher  latitudes  by 
the  increasing  correlation  between  neighbors  on  the  latitudeAongitude  grid,  are  more  definitive  in 
filtering  out  the  high-degree  components.  Special  weighting  functions  (such  as  the  Gaussian 
function)  can  further  reduce  the  aliasing  that  remains  due  to  the  ringing  of  the  uniformly  weighted 
average. 

These  results  were  presented  (Jekeli,  1994)  and  accepted  for  publication  (Jekeli,  1995)  with  the 
following  summary.  Presented  material  and  a  preprint  of  the  submitted  paper  are  included  in  the 
Appendix  to  this  report. 


5 


Summary 

The  methods  of  harmonic  analysis  on  the  sphere  were  investigated  from  the  point  of  view  of 
aliasing.  The  aliasing  error  was  formulated  in  terms  of  spherical  harmonic  coefficients  in  the  case 
when  the  function  is  sampled  on  a  regular  grid  of  latitudes  and  longitudes.  This  yielded  several 
results,  some  of  which  are  known,  but  perhaps  more  analytically  illustrated,  here:  1)  The  simple 
quadratures  method  and  related  methods  are  biased  even  with  band-limited  functions.  2)  A  new 
method  that  eliminates  this  bias  is  also  superior  to  Colombo's  method  of  least  squares  in  terms  of 
reducing  aliasing,  because  it  estimates  more  coefficients  at  little  added  computational  cost.  3)  But, 
a  simple  modification  of  the  model  makes  the  least  squares  method  identical  to  the  new  method  as 
one  is  the  dual  of  the  other.  The  new  method  solves  for  as  many  coefficients  as  data  and,  in  the 
presence  of  data  noise,  it  need  not  use  fabricated  noise  covariances  to  make  the  solution 
numerically  tractable.  4)  The  essential  elimination  of  aliasing  can  only  be  effected  with  spherical 
cap  averages,  not  with  the  often  used  constant  angular  block  averages. 


3.  Personnel 

Persons  who  carried  out  activities  under  this  contract  were:  Richard  H.  Rapp,  Christopher 
Jekeli,  Nagarajan  Balasubramania,  and  Dru  Smith.  Professor  Rapp  and  Jekeli  are  professors  in  the 
Department  of  Geodetic  Science  and  Surveying  while  Dr.  Balasubramania  and  Mr.  Smith  were 
Graduate  Research  Associates. 

4.  Travel 

During  this  contract  Professor  Rapp  attended  three  meetings  at  the  Goddard  Space  Flight  Center 
and  one  meeting  at  the  Defense  Mapping  Agency  Aerospace  Center.  Professor  Jekeli  attended  two 
meeting  at  the  Goddard  Space  Flight  Center. 

5.  Publications  and  Presentations 

The  following  list  summarizes  the  publications  and  presentations  that  were  made  describing 
activities  carried  out  under  this  project: 

Balasubramania,  N.,  Definition  and  Realization  of  a  Global  Vertical  Datum,  Phillips  Laboratory 
Scientific  Report  No.  1,  PL-TR-94-2144,  April  1994;  also  Report  No.  427,  Department  of 
Geodetic  Science  and  Surveying,  The  Ohio  State  University,  Columbus,  117p.,  April  1994, 
ADA283303. 

Jekeli,  C.,  Harmonic  analysis,  aliasing,  and  the  mean  anomaly.  Presented  at  the  meeting  of 
Working  Group  I,  Combination  Methods  for  High-Degree  Expansions,  NASA/DMA  Joint 
Gravity  Field  and  Geoid  Improvement  Project,  4  October  1994,  NASA,GSFC,  1994. 

Jekeli,  C.,  Spherical  harmonic  analysis,  aliasing,  and  filtering.  In  press,  Manuscripta  Geodaetica., 
1995 

Rapp,  R.H.,  A  World  Vertical  Datum  Proposal,  prepared  for  lAG  Special  Study  Group  5.149, 
lOS  Deacon  Laboratory,  Wormley,  England,  December  1993 

Rapp,  R.H.,  A  Global  Vertical  Reference  Frame,  prepared  for  Geodetic  Fixing  of  Tide  Gauge 
Bench  Marks,  meeting,  Wormley,  England,  December  1993 
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Rapp,  R.H.,  Vertical  Datums  in  Land  and  Marine  Positioning,  in  Prod,  of  International 
Symposium  on  Marine  Positioning,  Hannover,  p.  95-99,  Marine  Technology  Society,  1994 

Rapp,  R.H.  with  R.S.  Nerem,  A  Joint  GSFC/DMA  Project  for  Improving  the  Model  of  the 
Earth’s  Gravitational  Field,  presented  at  the  meeting  of  the  International  Geoid  Commission, 
Graz,  Austria,  September  1994 

Rapp,  R.H.,  The  Use  of  Potential  Coefficient  Models  in  Computing  Geoid  Undulations,  prepared, 
in  part,  for  the  International  School  for  the  Determination  and  Use  of  the  Geoid,  October,  1994; 
to  appear  in  the  Bulletin  Series  of  the  International  Geoid  Service,  1995 


Appendix 


This  appendix  contains  the  following  material  prepared  under  this  contract: 


Rapp,  R.H., 
Rapp,  R.H., 
Jekeli,  C., 
Jekeli,  C., 


A  World  Vertical  Datum  Proposal,  December  1993 
Vertical  Datums  in  Land  and  Marine  Positioning,  September  1994. 
Harmonic  Analysis,  Aliasing,  and  the  Mean  Anomaly,  October  1994. 
Spherical  Harmonic  Analysis,  Aliasing,  and  Filtering,  1995. 
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A  World  Vertical  Datum  Proposal 

Richard  H.  Rapp 
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December  17, 1993 


Aspects  of  Vertical  Datum 

•  Origin  of  Vertical  Datums  Associated  with 
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Several  Hundred  Vertical  Datums  in  the  World 


Vertical  Datums  in  Practice 

•  Country  or  Area  Datums 
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Scientific  Datum 


Vertical  Datum  Problems _ 

Regional  Datums  Have  Different  Reference  Surfaces 
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Is  Not  Available  in  All  Areas 


Regionat  Datum  Connections 

h  =  elliDsoidal  heiaht 
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Given  h  values  on  related  datums,  values  of  N,  B 
values  can  be  determined. 


_ Data  Used  for  Tests _ 

•  Approximately  2000  Doppler  positioned  stations 
from  1971-1986 
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(Units  are  cm) 


Selected  Connection  Results 

Germanv/Enqland  Datum  Difference 
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Ideal  World  Height  System  (WHS) 
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Difficulty  in  acquiring  well  distributed  and 
sufficiently  accurate  data 


A  Simplified  World  Height  System 

h  =  ellipsoidal  height  N  =  geoid  undulation 

H  =  orthometric  height  ^  =  height  anomaly 
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Given  h  and  ^  calculate  H 


What  is  the  Reference  Surface? 
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Accuracy 


Accuracy  of  Geoid  Determination  - 1995 
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World  Height  System  1995(?) 

This  WHS  will  be  based  on  satellite  positioning 
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•  Transformation  functions  can  be  developed  to 
convert  the  WHS95  system  to  local  datums  or  vice 
versa  when  needed. 
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representation  of  the  Earth's  gravitational  potential. 
Such  an  effort  is  now  being  planned  as  a  joint 
NASA/DMA  project. 


Vertical  Datums  in  Land  and  Marine  Positioning^ 

o 


Richard  H.  Rapp 

Department  of  Geodetic  Science  and  Surveying 
The  Ohio  State  University 
Columbus,  Ohio  43210 


ABSTRACT 


Vertical  datums  provide  the  reference  surface  for  measuring  heights  of  the  topography. 
The  heights  are  usually  expressed  as  mean  sea  level  heights  or  orthometric  heights  or  normal 
heights,  etc.  The  reference  surface  for  the  vertical  datum  is  ideally  an  equipotential  surface,  the 
geoid,  of  the  Earth's  gravity  field.  Since  sea  level  (or  mean  sea  level)  deviates  from  the  geoid  by  1 
to  2  meters,  the  various  vertical  datums  of  the  world  have  different  reference  surfaces  for  their 
origin.  This  presentation  will  first  consider  the  height  differences  between  the  reference  surfaces 
of  some  selected  (Australia,  United  States,  England,  Germany)  vertical  datums.  The  discussion 
will  then  turn  to  the  need  for  a  unique  vertical  datum  in  land  areas  and  then  ocean  regions. 
Whether  one  considers  land  or  ocean  applications,  a  case  will  be  made  that  the  geoid  should  be 
adopted  as  the  reference  surface  recognizing  that  the  absolute  portions  of  the  geoid  can  only  be 
determined  to  an  absolute  accuracy  of  10  to  100  cm  depending  on  data  available  for  use.  To 
achieve  such  accuracy  at  aU  wavelengths,  improvements  in  our  biowledge  of  the  Earth's  gravity 
field  are  needed.  Some  improvements  can  be  made  processing  existing  or  soon  to  be  avdlable 
data.  •  • 

Introduction 

Heights  and  depths  are  a  common  quantity  of  vast  usefulness.  Historically  height 
determination  was  referenced  to  mean  sea  level  which  was  accessible  through  tide  gauge 
measurements.  Heights  were  then  given  with  respect  to  mean  sea  level.  Although  one  finds 
normal  heights  and  orthometric  heights  in  use  today,  they  have  a  common  need  for  a  suitable 
reference  surface.  As  positioning  needs  evolved  numerous  vertical  damms  were  developed  for 
various  countries,  regions,  states  within  countries,  etc.  Today  there  are  100-200  different  vertical 
datums  in  the  world. 

The  concept  of  the  datum  is  to  have  a  gravity  equipotential  surface  from  which  orthometric 
heights  or  geopotential  numbers  can  be  measured.  Since  mean  sea  level  lies  on  different 
equipotential  surface  at  different  locations,  those  datums  based  on  mean  sea  level  are  not  based  on 
the  same  equipotential  surface.  Consequendy,  the  heights  from  different  datums  are  inconsistent  at 
the  ±2  meter  level. 

In  addition,  one  must  recognize  there  are  some  countries  that  do  not  have  access  to  an 
ocean  to  determine  a  mean  sea  level.  Such  areas  thus  must  bring  to  their  region  level  lines  from 
other  areas  having  such  access. 

Interest  in  smdying  the  relative  locations  of  vertical  datums  and  the  development  of  a  global 
vertical  datum  has  been  around  for  many  years  (Mather^  Rizos,  Morrison,  1978;  Colombo,  1980; 
Rapp,  1983;  Heck  and  Rummel,  1990).  More  accurate  space  positioning  and  improved 
knowledge  of  the  Earth’s  gravity  field  has  made  it  possible  to  study  the  differences  in  some  major 


^presented  at  the  International  Symposium  on  Marine  Positioning,  Hannover,  Germany,  September  1994 
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vertical  datums.  Two  recent  studies  are  those  of  Rapp  (1994)  and  Balasubramania  (1994).  The 
relationship  between  the  French  and  British  vertical  datums  has  been  discussed  by  Willis  et  al. 
(1989).  The  vertical  datums  in  the  Baltic  region  have  been  studied  by  Kakkuri  (1994)  and  Pan  and 
Sjoberg  (1994).  The  Rapp  smdy  used  Doppler  positioned  stations  while  the  Balasubramania  study 
used  stations  more  accurately  positioned  through  satellite  laser  measurements.  _  An  example  of 
results  is  shown  in  Figure  1  which  corresponds  to  Figure  Al  of  the  Balasubramania  report 


WHS  W  =  Wo 


Ref.  Ellipsoid  U  =  Uq 


Figure  1.  Vertical  Datum  Relationship  Based  on  Space  Located  Stations.  Units  are  cm. 

This  figure  shows  the  differences  between  selected  vertical  datums  and  the  ideal  reference 
surface.  Also  shown  is  the  reference  ellipsoid  (a=  6378136.3  m)  adopted  for  these  computations. 
One  sees,  for  example,  the  origin  difference  between  NAVD88  and  the  Australian  datums  (AHD) 
is  150  cm.  The  25  cm  difference  between  the  ideal  surface  and  the  adopted  reference  ellipsoid 
implies  that  a  better  radius  for  the  ellipsoid  would  have  been  6378136.55  m  which  is  consistent 
witii  what  is  being  found  with  Topex  altimeter  data. 

Although  much  attention  in  the  past  has  focused  on  vertical  damms  for  land  applications 
analogous  problems  occur  in  the  marine  environment  when  bathymetry  data  and  nautici  charts  are 
considered.  In  the  ocean  areas  there  are  also  numerous  datums  used  in  referencing  a  depth.  These 
damm^  are  usually  associated  with  some  specific  tidal  situation.  A  recent  review  is  given  is  Kumar 
(1994). 
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Global  Height  Systems 

As  the  precision  of  positioning  improves  with  the  use  of  GPS  it  becomes  imponant  to 
consider  the  determination  and  implementation  of  a  global  height  system.  Recent  studies  (e.g.  Xu 
and  Rummel,  1991;  Rapp  and  Balasubramania,  1992)  have  examined  ways  to  consider  the  global 
vertical  datum  problem.  Two  types  of  solutions  are  discussed.  One  combines  a  variety  of  space 
and  tenrestrial  data  to  define  an  optimum  system  connected  to  mean  sea  level  averaged  over  the 
oceans.  Such  data  would  include  geocentric  station  positions,  mean  sea  level  defined  at  tide  gauge 
stations,  gravity  field  information,  satellite  altimeter  data,  and  dynamic  ocean  topography. 
Although  concepmally  the  merger  of  such  data  is  desirable,  there  are  numerous  difficulties  in  terms 
of  data  coverage  and  data  accuracy  that  may  prevent  the  ideal  solution  from  being  realized. 

An  alternate  approach  is  to  use  the  geoid  as  the  reference  surface  and  let  this  geoid  be 
defined  by  a  set  of  geoid  undulations  derived  gravimetrically.  These  undulations  would  then  be 
given  with  respect  to  the  best  (as  of  some  epoch  and  agreement)  ellipsoid.  The  heights  (h)  above 
this  ellipsoid  would  be  determined  from  space  techniques,  primarily  GPS.  Knowing  h  and  N  the 
orthometric  height  can  be  determined: 

H  =  h-N  (1) 

If  normal  heights  (H*)  are  desired  the  gravimetric  information  could  be  used  to  calculate  the  height 
anomaly  (C)  so  that: 

H=^  =  h-C  ■  (2) 

In  these  computations  care  must  be  taken  to  be  in  the  proper  reference  system  or  frame.  This 
system  could  be  the  satellite  frame  defined  by  WGS84  and  the  constants  of  the  WGS84  ellipsoid. 
The  height  determination  implied  by  (1)  or  (2)  can  also  be  applied  to  height  difference 
determinations  knowing  that  Ah  can  be  very  accurately  determined  through  GPS  procedures. 

Similar  procedures  could  take  place  in  ocean  areas  for  bathymetry  mapping.  In  the  case  of 
the  ellipsoid  height,  suitably  corrected  for  tidal  effects,  the  ship  could  be  converted  to  a  height  with 
respect  to  the  geoid.  The  depths  measured  by  the  ship  could  then  be  referred  to  the  geoid  &umar, 
1994). 

Geoid  Undulation  Determinations 

The  usefulness  of  this  method  of  using  the  geoid  as  the  datum  origin  resides  on 
determining  geoid  undulations  to  a  sufficient  accuracy.  Today  geoid  undulations  can  be  calculated 
from  potential  coefficient  models  to  degree  360  (Rapp,  Wang,  Pavlis,  1991)  or  the  combination  of 
such  models  with  terrestrial  gravity  data.  The  accuracy  of  such  determinations  varies  from  ±30  cm 
in  ocean  areas  to  ±5  m  in  areas  lacking  terrestrial  gravity  data.  This  accuracy  could  be  improved  if 
better  satellite  and  terrestrial  data  is  used  in  the  development  of  degree  360  models.  Considering 
current  results  one  could  expect  geoid  undulations,  from  a  new  degree  360  to  achieve  an  accuracy 
from  ±25  cm  in  ocean  areas  to  ±1 .5  m  in  land  areas  not  adequately  covered  with  gravity  data. 

The  Joint  Project 

Recognizing  the  need  for  precise  geoid  undulation  determinations  for  both  geodetic  and 
oceanographic  applications,  the  Defense  Mapping  Agency  and  NASA/Goddard  Space  Flight 
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Center  have  started  a  project  to  determine  a  significantly  improved  Earth  gravity  model  based  on  a 
spherical  harmonic  representation  to  degree  360.  The  new  model  will  use  terrestrial  gravity  data 
from  areas  not  previously  available  and  from  revised  and  updated  files  in  regions  currently 
represented  in  existing  models.  New  satellite  tracking  data,  including  GPS  tracking,  will  be 
incorporated  in  the  new  model.  Satellite  altimeter  data  from  Topex/Poseidon,  Geosat,  and  ERS-1 
will  dso  be  used. 

Activities  on  the  development  of  the  new  model  are  now  underway.  Test  solutions  will  be 
started  in  October  1994  with  final  data  sets  expected  in  July  1995.  Candidate  final  models  should 
be  available  in  October  1995  with  selection  of  the  final  model  by  early  1996.  A  description  of  the 
techniques  and  data  to  be  used  is  given  by  Rapp  and  Nerem  (1994). 

The  output  of  the  joint  project  will  include  the  potential  coefficient  model  and  a  gridded 
geoid  undulation  files  based  only  on  the  degree  360  model.  These  gridded  undulations  can  be  used 
for  the  determination  of  orthometric  heights  (or  height  differences)  from  ellipsoidal  heights  derived 
from  GPS.  Care  will  be  taken  to  have  the  geoid  undulations  in  the  WGS84  system. 

Conclusions 

A  practical  implementation  of  a  world  height  system  or  vertical  datum  is  being  considered 
using  the  geoid  as  the  conceptual  reference  surface.  Although  this  surface  will  not  be  directly  tied 
to  mean  sea  level,  it  will  be  associated  with  a  global  mean  sea  level  through  the  proper  choice  of  the 
parameters  of  the  ideal  ellipsoid.  The  key  element  in  the  concept  will  be  in  determining  geoid 
undulations  to  meet  users  needs.  Initially  this  will  be  done  by  developing  a  new  degree  360 
potential  coefficient  model  based  on  significantly  new  data  sets.  This  will  enable  undulations  to  be 
derived  for  most  regions  to  an  accuracy  of  ±1.5  m  or.  better.  Improved  determinations  can 
subsequently  take  place  for  specific  regions  through  detailed  geoid  undulations  computations 
incorporating  high  resolution  terrestrial  gravity  data. 

The  geoid  undulations  will  play  a  key  role  in  converting  ellipsoidal  heights  into  orthomeiric 
heights  with  respect  to  the  single,  world  wide  defined,  reference  surface,  the  geoid.  This  will  be 
important  for  both  land  applications  and  marine  charting  purposes.  For  the  first  time,  one  has  the 
prospect  of  having  a  single  vertical  datum.  The  use  of  this  single  datum  by  all  countries  would  not 
be  realistic.  However,  the  development  of  transformations  between  a  local  datum  and  the  global 
datum  is  quite  possible.  This  would  enable  the  conversion  of  heights  from  a  local  vertical  datum  to 
the  global  datum  or  vice  versa. 

The  increased  accuracy  of  positioning,  both  on  land  and  in  the  oceans,  makes  the 
development  of  a  world  height  system  a  natural  progression  from  a  determination  of  a  world 
geodetic  system  (e.g.  WGS84).  Steps  are  now  being  taken  to  have  information  to  implement  the 
world  height  system/vertical  datum  concept  in  1996. 
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HARMONIC  ANALYSIS,  ALIASING,  AND 
THE  MEAN  ANOMALY 


O) 

D) 

C 

’0 

_0 

< 

• 

c 

0 

‘0 

> 

0 

E 

0 

0 

> 

c 

0 

< 

'0 

>. 

-“ti 

0 

> 

0 

C 

v_ 

< 

0 

,o 

C 

’c 

0 

o 

0 

E 

0 

0 

X  . 

O 

0 

0 

0 

0 

0 

*3 

■G 

O 

0 

• 

• 

o 

CD 


O 

X) 

E 


^  •- 


CO  ^ 


.p 


£ 

p 

-c 

§- 

p 

-c 


0 

0 

0 

x: 

Q. 

O 

4— » 
(/) 


o 


CO 

c 

■>s 

(D 

CD 

“O  >. 
C  .tii 
03  CO 

CD  CD 
O  > 

.1  3 

o 

CD  ^ 
o  ^ 

■•s  ^ 

•O  o 

^  *r" 

^  o 

CD 


ZD 

TD 

0) 

O* 

O  B 

T—  « 

CVI 

c?o 


CO 

o 

cc 

CO* 

13 

c 

D> 

cc 

E 


o 


e 

0 

E 

tr 

03 

Q- 

0 

Q 


X 

o 

0 
=3 
Xi 

E  = 

03 

£ 


0 

jD 

O 

H— » 

O 

O 


0 

0 

.3 

0 

0 


■  • 

>> 
■4— > 

3 

O 


Q 


29 


Spherical  Harmonic  Analysis  of  Gravity  Anomalies 
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Eliminate"  or  Reduce  Aliasing  by  Averaging  Gravity  Anomaly 
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Models  of  Harmonic  Analysis  in  Practice 
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Ag(0,>\.)  =  X  X  Tn,m  Yn,m(0,^)  OiTiits  Higher-Degree  Spectrum 

m=-N  n=jm| 


Residual  Geoid  Signal  According  to  Kauia's  Rule 

[64/N,  meters] 

N  =  360 


Quadratures  Harmonic  Analysis 


degree 


Conventional  Least-Squares  Harmonic  Analysis 


degree 


Options  to  Model  the  Mean  Anomaly 
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Gaussian  vs.  Pellinen  Smoothing  Factors, 
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latitude  [degrees] 


•  Simple  quadratures  (inch  variations)  are  deficient,  in  theory,  for  harmonic  analysis 
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-  Gaussian  filter  is  better  than  Pellinen  filter  (unweighted  average) 

-  least-squares  collocation  is  not  necessarily  the  best 


Spherical  Harmonic  Analysis,  Aliasing,  and  Filtering 
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Abstract.  The  currently  practiced  methods  of  harmonic 
analysis  on  the  sphere  are  studied  with  respect  to  aliasing 
and  filtering.  It  is  assumed  that  a  function  is  sampled  on 
a  regular  grid  of  latitudes  and  longitudes.  Then, 
transformations  to  and  from  the  Cartesian  plane  yield 
formulations  of  the  aliasing  error  in  terms  of  spherical 
harmonic  coefficients.  The  following  results  are  obtained: 
1)  The  simple  quadratures  method  and  related  methods  are 
biased  even  with  band-limited  functions.  2)  A  new 
method  that  eliminates  this  bias  is  superior  to  Colombo’s 
method  of  least  squares  in  terms  of  reducing  aliasing.  3) 
But,  a  simple  modification  of  the  least-squares  model, 
makes  it  identical  to  the  new  method  as  one  is  the  dual  of 
the  other.  4)  The  essential  elimination  of  aliasing  can 
only  be  effected  with  spherical  cap  averages,  not  wiUi  the 
often  used  constant  angular  block  averages. 


Introduction 

The  harmonic  (or  spectral)  representation  of  functions 
on  the  sphere  has  proved  useful  in  modeling  the  gravity 
field,  the  magnetic  field,  the  topography,  the  ocean  tides, 
etc,  not  just  for  the  Earth,  but  for  any  roughly  spherical 
planet  (Rapp  1977;  Schmitz  and  Cain  1983;  Kaula  1993). 
Such  harmonic  representations  enable  a  more  precise 
interpretation  by  wavelength  of  their  magnitudes,  their 
coherence,  and  their  measurability.  Usually,  these 
functions  are  richly  endowed  with  harmonics  to  a  high 
degree,  decaying  relatively  slowly  with  increasing  detail. 
For  example,  the  Earth’s  gravity  field,  after  the  first  few 
harmonics  (which  certainly  dominate),  decays 
approximately  linearly  on  the  logarithmic  scale,  and 


today's  representations  include  over  130  thousand  terms 
(Rapp  and  Pavlis  1990). 

It  is  a  significant  consequence  of  this  fine  structure  in 
the  functions  that  the  harmonic  components  detennined 
from  a  finite  set  of  function  values  are  biased  estimates  of 
the  true  components.  This  is  a  well  known  phenomenon 
in  signal  analysis,  where  even  though  one  applies 
orthogonal -operators  to  the  data,  they  cannot  filter  out 
harmonics  finer  in  detail  than  that  dictated  by  the 
sampling  interval.  This  is  known  as  "aliasing.”  Another 
interpretation,  though  less  precise,  is  that  of  attempting  to 
determine  all  (possibly  an  infinite  number  of)  harmonics 
from  a  given  number  of  function  values.  In  the  general 
case,  this  is  an  underdeiermined  problem  and  the  solution 
invariably  lumps  together  harmonics  in  a  prescribed  way. 
Spectral  aliasing  in  functions  defined  on  the  line  or  on  the 
Cartesian  plane  is  well  understood;  it  is  less  obvious  on  a 
non-Euclidean  surface  such  as  the  sphere. 

The  spherical  harmonic  analysis  is  examined  here 
with  the  object  of  clearly  understanding  the  effect  of 
aliasing  in  conventional  techniques,  as  well  as  in  the 
modem  techniques  developed  by  Colombo  (1981)  that  are 
used  to  compute  today’s  models.  In  so  doing,  a  new 
formula  for  the  analysis  is  rendered  that  reduces  the 
aliasing  effect  better  than  the  least- squares  approach.  It  is, 
however,  also  shown  that  the  new  technique  is  equivalent 
to  a  modification  of  the  least-squares  model,  requiring 
some  extra  computations.  The  results  also  indicate  clearly 
which  types  of  averaging  (filtering)  do  and  do  not 
eliminate  aliasing.  It  is  assumed  throughtout  that  the  data 
are  provided  discretely  at  every  node  of  a  uniform  grid. 
(Aliasing  in  the  case  of  a  nonunifonn  data  distribution 
was  studied  by  Sans6  1990.)  Though  the  errors  in 
measuring  the  function  on  the  sphere  are  not  considered 
initially,  their  effect  on  the  derived  procedures  is  discussed 
in  conclusion. 
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Preparatory  Concepts 

A  periodic  function  iniegrable  over  its  period  may  be 
represented  either  in  terms  of  its  independent  (space  or 
time)  variable  or,  where  it  is  continuous,  in  terms  of  its 
Fourier  spectrum,  being  the  set  of  coefficients  in  its 
corresponding  Fourier  series  (Priestley  1981): 


g(x)=  X  Gkc’-’''^  ,  xe;^ 

k=>*o 

(la) 

(•T 

Gi  =  —  gW  dx  ,  ke 

(lb) 

Jo 


where  r  =  -1,  T  is  the  period  (g(x+T)  =  g(x)),  !^is  the  set 
of  real  numbers,  Z  is  the  set  of  integers,  and  the  integer  k 
may  be  termed  ’’frequency,”  or  "wavenumber."  In  practice, 
for  example  from  measurements,  one  knows  only  discrete 
values  of  g  and  the  corresponding  spectrum  is  the  Disaete 
Fourier  Transform  (DFT): 


The  infinite  sum  on  the  far  right  in  (4)  is  an  alias  of 
for  any  frequency,  k,  such  that  Ikl  <  kj^  =  K/2,  where  is 
the  so-called  Nyquist  frequency.  The  deteimination  of  the 
spectrum  {G^}  of  the  function  g  from  the  disaete  sequence 
{g^},  using  (2b),  is  subject  to  an  aliasing  error  as 
formulated  in  (4).  Note  that  if  the  function  g  does  not 
contain  harmonics  with  frequency  above  the  Nyquist 
frequency,  then  there  is  no  aliasing  error.  Such  functions 
will  be  called  band-limited.  It  is  further  noted  that  a 
function  with  infinite  bandwidth  can  be  filtered  to  make  it 
band-limited.  This  procedure  is  discussed  later  with 
respect  to  the  particular  application  of  spherical  harmonic 
analysis. 

Obviously,  similar  formulas  hold  for  periodic 
functions  defined  on  the  Cartesian  plane.  They  are 
specialized  here  to  the  case  where  g  is  a  real  function  of 
two  variables,  0  and  X,  and  is  periodic  in  0  with  period  n 
and  in  X  with  p^od  2jc.  Then 


gs  =  X  Gk  ,  s  =  0 . K-1 

k=0 

(2a) 

g(9,X)  =  X  X  Gtm  e,Xe  3^ 

k=-^  n)=>^ 

(5a) 

o  ,  K-1 

Gem - 

2:t'J 

f’-' 

Gi  =  A2L  y  gj  6--’''=^  .  k  =  0 . K-1 

T  sTo 

(2b) 

[  J 

g(e,X)e-'^''®*'^>d6dX  . 

'o 

k,me  Z 

where  it  is  assumed  that  the  function  \'alues  are  known  on 


a  regular  grid,  i.e.,  one  with  constant  spacing.  Ax  =  T/K. 
The  unsymmetric  transform  (2a,b),  where  k  is 
nonnegative,  is  preferred  here  as  being  more  conventional 
and  avoiding  the  separate  treatment  of  even  and  odd  K. 
Since  the  DFT  is  periodic  with  frequency  period  K: 

Gk+K  =  Gic ,  for  any  k  (2c) 


where,  with  *  denoting  the  complex  conjugate, 

(jTc.-m  ==  G-k.m  (^) 

For  a  regular  grid  with  constant  spacings  given  by  A0  = 
n/K  and  AX  =  2kM,  where  0^  =  sA0,  X^  =  tAX,  the 
Discrete  Fourier  Transform  pair  is 


it  is  noted  that  equations  (2a,b)  are  equivalent  to  the 
symmetric  forms,  where  k  (and  s)  ranges  from  -(K-l)/2  to 
(K-l)/2  (if  K  is  odd).  In  the  notation  of  (2a,b),  the 
frequencies  from  (K+l)/2  to  K-1  are  used  to  represent  the 
negative  frequencies  from  -(K-l)/2  to  -1. 

The  relationship  between  the  spectra  (lb)  and  (2b)  can 
be  found  by  substituting  x  =  sAx  into  (la): 

g5  s  g(sAx)  =  X  Qi  e*^*'*^ 

7i  - 

=  I  i 

k=0  j=— 

from  which,  by  comparing  with  (2a),  one  has  for  any  k, 

Gk  =  X  Q'K+k  =  Gk  +  X  GjK+k  (4) 


X-l  M-l  o 

gs.t  =  gOs-W  =  £  X  Gun  . 

k=0  m=0 

(6a) 

s  =  0,...,K-1  ,  t  = 


O  ^  K-1  M-l 

Gua  =  -LX  Xgs.e-*W*tnVM)  ^ 

KM  s=o  1=0 

(6b) 

k  =  0,...,K-1  ,  m  = 


Because  the  values  g^  are  real, 

O  0  ^ 

GkJ^m  =  GK-kjn  (6c) 

which  explains  the  limit  of  interest  on  m.  M  is  assumed 
to  be  even,  wthout  Joss  in  generality.^  By  the  peri^city 
in  k  and  m.  Go jn  =  GoaM-m  and  Gk,o  =  Gx-to,  while  Go.o  is 
real.  Hence,  the  complex,  discrete  spectrum,  (Gkjn), 
contains  only  KM  independent  coefficients.  Analogous  to 
(4),  the  relationship  between  the  total  spectrum  (5b)  and 
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the  sample  spectrum  (6b)  is,  for  any  k  and  m, 

O  OO  oo 

Gkjn  =  Glc^m  +  X  X  GuK4t.vM+m  (7) 

The  double  sum  in  (7)  represents  the  aliasing  error  for 
frequencies  iki  <  KJ2  and  Im!  <  M/2. 

Instead  of  the  plane,  the  real  function  g(0,X)  can  also 
be  defined  on  a  unit  sphere  and  then  expressed  in  lenns  of 
a  series  of  spherical  harmonics  rather  than  sinusoids. 
There  are  many  variations  in  such  a  representation.  For 
simplicity,  the  complex  exponential  form  is  used  here  and 
its  relationship  to  the  more  customary  form  is  also  given 
for  the  reader's  convenience.  The  spherical  harmonic 
series  defines  the  complex  Legendre  spectrum  {y  }  of 
the  function  as  follows: 

g(6,X)=X  X  Yfl.m  Pn.|ra^COS9)  e*“^  ; 

(8a) 

0<e<7C,  0<X<27r 


to  the  coefficients  of  Heiskanen  and  Moritz  (1967,  p.31, 
(T76))by 

i ~  i  bn jn)  >  m  >  0 
yQ.m  =  /aD,o,  m  =  0  (10) 

[^iQ,m  +  i  bnjn) ,  m  <  0 

Because  of  (8d),  one  need  only  solve  for  coefficients  with 
non-negative  orders  although,  formally,  the  series  must 
include  m  <  0. 

The  Fourier  Spectra  of  _  and  sin0  P_  _ 

^  n,m  n,m 

It  will  be  necessary  to  know  the  Fourier  spectra  of 

and  of  sin0-  both  defined  for  0  <  9  < as 

usual,  and  for  0  <  0  and  0  >  tc  by  periodic  extension. 

This  means  that  the  number  of  frequencies  of  the  spectrum 

in  each  case  is  infinite^  even  though  these  functions  are 

polynomials  of  sinusoids  (Hobson  1965,  p.96).  For  0  < 

0  <  Tt,  and  any  integer,  k,  let 

PDjn(0)  =  pD.m(cOS0)  and  Pa.m(0  +  kjx)  =  Pn,in(0) 


Yn.m 


4;iem 


Jo 


g(8,X)  Pq.:i:3.(cos9)  sin0  d0  dX  ; 

.0 


(11a) 


SDjn(0)  =  Sin9  Pajn(COS0)  and  Sn.m(9  k>’t)  =  Sn.mO) 


where 


n>0,  -n<m<n 


(lib) 


(8b) 


In  accord  with  (la)  and  (lb),  the  corresponding  Fourier 
transform  pairs  are 


£m  — 


m  =  0 
m^O 


(8c) 


Pa^(e)  =  ^X  ..  06 


(12a) 


and,  since  g  is  real. 


Yn.-m  — 


m 


The  integers  n  and  m  are,  respectively,  the  degree  and  order 
of  the  harmonic  coefficient,  y  .  The  P  in  (8a)  an  {8b) 
are  caOed  associated  Legendre  functions  of  the  first  kind, 
defined,  for  the  present  purpose,  for  non-negative  integers 
n,m  with  m  <  n  and  for  0  <  0  <  k.  Furthermore,  they  are 
assumed  to  be  normalized  such  that 


1 

2 


rrc 

(PQjn(COS8))^  SUl9  d9  =  Em 

.0 


(9) 


(again,  for  simplicity,  the  conventional  over-bar  notation 
is  omitted).  The  complex  harmonic  coefficients  are  related 


=  l  P„j„(e)e-‘-''®d8  ,  k€Z  (12b) 

JO 

and 

SDjn(9)=  X  ,  0ei?^  (13a) 

=  sine  PojnO)  de  ,  ics  2  (13b) 

Jo 

It  is  noted  that  because  P  and  S  are  real  functions, 

ii,jn  Djn 

pr  =  (p^)*  and  ar  =  (a“f )‘  (14) 

The  orthogonality  of  the  associated  Legendre  functions 
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implies  that  the  infinite  sequences  {pt"}  and  [di]  are  Substituung  (12a)  and  (8a)  into  (5b)  and  taking  note  of 
orthogonal  for  n  j:  (18a,b),  one  obtains 


y  p*i'"o^'"  =  ?^6(j,n)  (15) 

where  5(j,n)  is  the  delta  function,  equal  to  zero  unless  j  = 
n,  when  it  is  one. 

Once  these  coefficients  are  computed  for  n,m  =  0,1 
and  for  ail  integers,  k,  of  interest,  the  remaining  spectra 
(n,m  >  1)  for  these  k  can  be  determined  using  well 
known  recursion  formulas  for  the  associated  Legendre 
function  (Abramowitz  and  Siegun  1972). 

The  discrete  inverse  Fourier  transforms  of  these 


functions  sampled  on  a  uniform 
{05  =  s  A0  :  s  =  0, K-1;  A0  =  n/K}  are 

grid 

k=0 

(16a) 

K-1 

C  -  V  J2:csk/K 

Onjn-.s  =  2-r  ^ 

k=0 

(16b) 

with  (see  (4)): 

Pk  =  2^  PuK+k  » 
u=— 

(17) 

The  discrete  spectra  in  (17)  are  also  amenable  to 
computation  by  recursion  formulas  (see  Appendix)  since 
the  left-hand  sides  are  simply  finite  sums  of  associated 
Legendre  functions.  Other  methods  to  compute  these 
include  simply  determining  the  Discrete  Fourier 
Transform  (A.l)  of  the  Legendre  functions  determined 
from  recursion  formulas,  or  substituting  their  equivalent 
finite  polynomial  of  sinusoids  into  (A.1). 

The  Relationships  Between  the  Fourier  and 
Legendre  Spectra 

The  following  orthogonality  relaiicxiships  hold: 


Gkjn  =  Yj  pt 

where  it  is  further  noted  that  y  =  0  for  1ml  >  n. 
Similarly,  substituting  (13a)  and  (5a)  into  (8b)  yields  Lee 
inverse  reladonship 

=  Gu^or'  (20) 

2Zm 

These  relationships  are  possible  because  of  the  definitions 
(lla,b).  The  alternative  suggested  by  Ricardi  and  Burrows 
(1972)  implicitly  assumes  periodicity  of  the  associated 
Legendre  functions  over  (0,2:1)  and,  though  their  spectral 
domains  are  then  finite,  in  this  form  they  cannot  be  used 
with  (5a)  and  (5b)  that  embody  periodicity  in  latitude  over 
(0,k).  On  the  other  hand,  Goldstein  (1978)  devised  an 
extension  of  the  function  g  onto  the  larger  domain,  0  <  9 
<  iTZy  0<X<  27C,  thus  taking  advantage  of  the  associated 
Legendre  functions  as  finite  polynomials  of  sinusoids.  In 
the  present  case  the  infinite  nature  of  the  their  spectral 
domains  will  not  cause  computational  problems  because 
in  the  end  only  the  discrete  spectra  (17)  are  needed. 
Albertella  and  Sacerdote  (1995)‘also  developed  procedures 
simOar  to  the  present  ones. 

Quadrature  Formulas 

Because  the  function  (e.g.,  the  gravin'  anomaly)  is 
supposed  to  be  given  on  a  regular  grid,  all  estimates  of  its 
Legendre  spectrum  have  the  form  of  so-called  quadrature 
fonnulas.  A  subset  of  these  is  the  discretization  of  the 
integral  (8b).  Numerous  variations  have  been  proposed; 
for  example,  the  simplest  is 

=  _iL_  I;  X  2s.t  PnanCcosGs)  sinGs  e^^^  (21) 

2£mKM  ^0  ^ 


2tz 


JL 

2n: 


^0 


if  m  m' 
if  m  =  m' 


and 


1 


gi20c*-k)8  ^0 


0,  ifk^k 

T  ,  if  k  =  k 


and  for  later  use: 


1=0  ■  1m  ,  j  =  m 


More  common  is  an  attempt  to  incorporate  the  fact  that 
(18a)  usually  the  data  are  averages  over  angular  cells.  In  this 
case,  the  Legendre  spectrum  of  the  mean  anomaly  is 
sought  and  the  spherical  harmonic  function  \-aiues  in  (21) 
are  replaced  by  their  averages  over  the  cells.  Either  way, 
(20)  or  its  alternative  with  averaged  spherical  harmonics  is 
only  an  approximation  and  leads  to  estimates  that  are 
(18b)  biased  due  to  aliasing  not  just  from  frequency  components 
higher  than  the  Nyquisi  limit,  but  from  components  of  ail 
degrees  greater  than  or  equal  to  the  order.  Indeed,  it  is  well 
known  (Pavlis  1988;  Albertella  et  al,  1993)  that  a  band- 
limited  function  analyzed  according  to  (21)  is  not 
(18c)  recovered  exactly  by  the  estimated  spectrum.  This  feature 
is  demonstrated  analytically  below,  but  only  for  the  case 
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exemplified  by  (21).  The  analogy  for  the  case  of  cell- 
averages  will  be  evident  because  the  integrals  of  the 
Legendre  functions  have  spectra  like  (17),  but  multiplied 
by  (sinkA0)/(kA0). 

With  the  aid  of  (13a)  and  (17),  equation  (21)  becomes 

=  (22) 

2£ni 

The  error  in  this  estimate  is  derived  by  substituting  (7) 
into  (22),  adding  and  subtracting  and  making  use  of 
(17),  (19)  and  (20).  Some  manipulations  yield 


2£:ni 


(23) 


K-l 


ZV*  V  °  Qjm| 

2^  2^  Pi:  fi^vM+a 


it=0  v=-**  j=:^vM+ml 

Now  shift  the  v=0  term  to  the  first  sum  to  get,  with  (15), 


The  first  sum  contains  exactly  K  harmonic  coefficients  for 
any  given  m,  and  these  are  to  be  solved,  given  One 
could  have  separated  the  first  two  sums  in  (26)  at  some 
integer  larger  than  Iml  +  K,  but  then  the  problem  is 
underdetermined  because  there  are  oniv  K  Fourier 

o 

coefficient  for  any  given  m. 

Let  Gm  and  y^  be^finite  K-vectors  containing, 
respectively,  coefficients  Giun  and  Yj  as  follows: 


yin  —[yirn  jn,  Y  ui  — »  yimj+K-ljn]  (2S) 


Also,  define  the  KxK  matrix 

Rm  — 

( 

... 

1 

*•  FK-1  J 

(29) 


yo.m  —  Yn.m  + 

i  (l  oT'  -  ^ 5(j.n))  -  (24) 

2£m  U=irT5|  \lc=0  ^  / 

i  i 

v=-«*  vM+m|  1:=0 

Clearly,  for  band-limited  data,  where  y^  ^  =  0,  n  >  M/2, 
the  last  set  of  sums,  where  j  >  M/2,  vanishes.  But,  the 
estimate  is  still  biased  ("aliased")  by  the  second  set  that 
contains  all  harmonic  coefficients  with  degrees  j  >  iml, 
including  y^^  The  analysis  procedure  derived  next 
shifts  the  total  error  of  the  estimate  to  frequencies  beyond 
the  Nyquist  frequency. 

A  New  Method 

Combining  (7)  and  (19),  and  using  (17)  leads  to 
Gta=  i:  i  ^k^''"'*'"'Yi.vM..  (25) 

v=:—  j=^v.M+ml 


This  is  the  true  model  linking  the  two-dimensional 
discrete  Fourier  transform  of  the  data  (6b)  to  the  actual 
Legendre  spectrum.  It  is  equivalent  to  the  model  (8a)  for 
the  nontransformed  data.  Again,  extracting  the  v  =  0 
contribution,  this  can  be  broken  into  three  parts: 


lm|  +-K-1  — 

G^=  I  ^rYi.m+  I  ^rYi.n 

j=iml  j=|mKK 


vi»0 


o^vM+mj 

Ft 


yj,vM+m 


(26) 


Then,  the  collection  of  equations  (26)  for  a  fixed  m  and  for 
k  =  0, K-l  is 

O 

Gm  “  Rm  ym  +  5Rm  5ym  +  ARm  A  /m  (-^0) 

where  5R  ,  5y  ,  AR  ,  and  Ay  are  semi-infinite  matrices 

m  *m  m  *m 

and  vectors  whose  elements  can  be  inferred  from  (26).  In 
particular,  5y^  contains  harmonic  coefficients  of  order  m 
and  desree  greater  than  Iml  +  K  -  1;  while  Ay^  contains 
all  coefficients  with  order  greater  than  or  equal  to  M/2. 
The  relative  domains  for  these  coefficients  are  shown  in 
Figure  1,  where  it  is  noted  that  K  is  completely 
independent  ofM. 


Fig.l;  Domain  of  spherical  harmonic  coefficients  estimated  by  (31). 
Aliasing  is  due  to  coefficients  Ay  and  Sy  in  their  respective  domains, 
see  eq.  (32). 
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Solving  (30)  for  yields  the  following  estimate  of 
the  harmonic  coefficients 

?n,  =  R;^Gn.  (31) 

with  an  aliasing  error  given  by 

Ym  =  “Rm  (SRm  ^Ym  +  ARm  AYm)  (32.) 

R  is  nonsinguJar  as  shown  later  (equation  (49)).  For  a 

IQ 

band-limited  function  the  aliasing  error  (32)  should 
vanish.  Band-limited  now  refers  to  two  Nyquist 
frequencies,  just  like  in  the  planar  case:  N(m)  =  M/2  and 
N(n)  =  m  +  K.  Thus,  the  spherical  analysis  of  functions 
having  harmonics  (n,m)  that  are  nonzero  only  for  n  < 
N(n)  and  m  <  N(m)  will  be  perfect,  without  aliasing, 
according  to  (31).  The  use  of  the  quadratures  formula  (21) 
would  introduce  errors  as  given  by  the  first  double  sum  in 
(24). 

It  is  convenient,  especially  in  efforts  to  eliminate  the 
aliasing  effect,  to  refer  the  Nyquist  limit  to  a  single 
frequency  in  terms  of  the  degree,  n.  This  is  N  = 
min(K,M/2).  where  usually,  of  course,  K=M/2,  This, 
however,  understates  the  number  of  recoverable  spectral 
componen^.  That  is,  from  the  KM  data,  g^  .,  that 
determine  Gm,  (31)  yields  exactly  KM  real  coefficients 
contained  in  the  complex  numbers  Ym.  0  <  m  <  M/2;  see 
Figure  1.  Note  that  for  each  such  m,  the  estimate  (31) 
requires  the  inversion  of  a  KxK  matrix. 


true  model,  only  an  approximation,  due  to  the  truncation. 
For  a  more  compact  notation,  define  the  following  K}vl- 
vector 

g  =  {gs.t} .  (35) 

the  (N  -  lml)-vector  (which  differs  from  the  K-vecior, 
by  the  number  of  elements) 

Ym  —  [Yni|,m?  Yni|+l,mj  •••>  YN-l.mj  »  (^^ 

and  the  KMx(N-Im[)  matrix 

Y.  .  { OT 

The  tildes  (-)  are  included  now  to  differentiate  between 
this  and  the  following  modification  to  the  approach.  With 
(35),  (36),  and  (37).  the  model  (34)  becomes 

m=-N+l 

Note  that  there  are  KM  data  (observations),  but  only 
approximately  N"  coefficients  (if  K=M/2,  then  N"  = 
KM/2),  The  true  model  (33)  can  be  written  like  (38)  as 
follows: 

g  =  "x  Ym  +  5Ym  S'/m)  Y.  ^ 


The  Least-Squares  Solution 

It  is  now  shown  that  the  procedure  above  is  similar  to 
the  least-squares  approach  developed  by  Colombo  (1981) 
and  implemented  by  Hwang  (1993),  However,  there  is  a 
subtle  difference  that  is  easily  eliminated.  Equation  (8a), 
with  summations  reversed,  is  rewritten  as  follows  to 
separate  symbolically  the  coefficients  to  be  estimated: 


g(es,xo  = 

I  I  +I  i  +  I  i  l(Yn.=>P^m(COSe,)e^) 

D=^m|  m^N+l  n=N  a=im|  j 

(33) 


The  usual  model  for  the  least-squares  approach  is  the 
truncated  spherical  harmonic  series,  being  the  first  part  in 
(33): 


N-l  N-1  .  . 

g$.t  =  z  X  Tiijn  Ptt.im.(cosej)  e“^  ; 

(34) 

0<s<K-l,  0<  t<M-l 


The  symbol  ~  is  used  to  emphasize  that  this  is  not  the 


where  5Ym  is  a  semi-infmite  matrix  with  KM  rows  and  its 
elements  are  spherical  harmonic  functions  with  degrees  n 
>  N.  5Ym  contains  coefficients  with  degree  n  >  N,  is 
the  same  as  in  (30)  and  AY  can  be  inferred  from  (33). 

By  the  discrete  orthogonality  of  the  complex 
exponential  (18c),  the  matrices  Ym  are  also  orthogonal: 

K7Yj  =  0,  m 

Multiplying  (38)  by  (Ym}  and  denoting  the  ncsmai  matrix 
by  Am  =  { Ym)  Ym  yields  the  "least-squares"  solution: 

Yn.  =  AmiYm)  g  (41) 

The  aliasing  error  is  obtained  by  substituting  the  true 
model  (39)  and  duly  noting  the  onhogonalities: 

^  =  Ym  +  (SYm5ym  +  AYm  AYa)  (42) 

Rgure  2  shows  the  parts  of  the  spectrum.  Sya  and  AYjj,. 
that  alias  the  estimated  part  (K=M/2  is  assumed).  Note 
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(46) 


the  similarities  of  the  equation  pairs  (41),  (42)  and  (31), 
(32);  however,  the  solution  (41)Js  aliased  more  than  the 
solution  (31).  Because  the  Am  are  (N4m!)x(N-lmI) 
matrices,  there  are  fewer  computations  associated  with 
(41)  (matrix  inverses  of  diminishing  dimension  as  m 
N-1). 


Substituting  (45)  into  (31)  and  equating  with  (44)  yields 
=  (47) 


Fig.2:  Domain  of  spherical  harmonic  coefficients  estimated  by  (41), 
the  conventional  least-squares  approach.  Aliasing  is  due  to 
coefficients  Ay  and  5y  in  their  respective  domains,  see  cq.  (42), 


A  Modification  To  The  Least-Squares  Model 

By  now  it  is  evident  that  the  difference  between  (41) 
and  (31)  lies  in  the  model.  Augmenting  (34)  to  include 
degrees  up  to  1ml  +  K  - 1,  the  model  becomes 

M/2 

g=  X  YmYm  (43) 

m=-M/2+l 

where  is  similar  to  Ym  but  now  each  has  K  columns 
instead  of  (N  ';;_Imi),  and  is  given  by  (28). 
Correspondingly,  Am  becomes  the  KxK  matrix  A^;  and 
the  least-squares  solution  (41)  becomes 

Ym  =  Am  (Yih)  g  (44) 

Since  as  many  coefficients  are  computed  as  data  given, 
this  is  no  longer  a  "least-squares"  solution.  If  g  is  band- 
limited,  then  both  (31)  and  (44)  yield  the  true  spectrum, 
which  means  that  (31)  and  (44)  are  identical  for  any  g. 
Hence,  in  general,  the  aliasing  error  in  this  modified 
approach  is  also  given  by  (32). 

_  o 

The  2-D  Fourier  transform  Gm  (K-vector)  is  a  linear 
operator  on  g: 

=  (43) 

where  Wm  is  a  KMxK  matrix: 


which  holds  for  arbitrary  g.  Therefore,  since 
Wm  (Wm7  =  KM  I,  where  I  is  the  identity  matrix, 

Am  =  (Y;f{w;fR™  (48) 

Finally,  it  can  easily  be  shown  that 
(Y;f{w;7=KMR;  ,  SO  that  (48)  becomes 

Am  =  KM  Rm  Rjn  (49) 

which  is  Paneval's  theorem  in  view  of  (29)  and  because 
the  elements  of  A  are 

m 

K-1 

^jjc  ^  Pjim{(C0S9s)  Pic Jir:((^Os6s)  (OO) 

s=0 

One  can  stan  with  (49),  being  Parseval’s  Tceorem, 
and  show,  first,  that  since  Am  is  a  full  rank  matrix,  Rm  is 
non-singular;  and,  second,  that  solutions  (31)  and  (44)  are 
equivalent.  Thus,  there  is  a  perfect  dualit>'  between  these 
solutions:  (44)  operates  directly  ^on  the  signal,  g;  (31) 
operates  on  its  Fourier  transform  G;  both  yield  identical 
results. 

Harmonic  Analysis  of  Averages 

In  order  to  make  use  of  all  gravity  data  around  the 
world  and  to  keep  the  computation  of  spherical  harmonic 
models  tractable,  usual  practice  is  to  transform  the  data 
into  mean  quantities,  i.e.,  averages  over  certain  defined 
blocks  or  areas  on  the  Earth’s  surface,  often  delineated  by 
constant  lines  of  latitude  and  longitude,  e.g.,  l®xl®  or 
30’x30'  blocks.  This  yields  a  new  function,  g,  in 
principle,  defined  for  every  point  on  the  sphere,  but 
sampled,  like  g,  on  a  regular  grid,  g  has  its  own  unique 
Legendre  spectrum,  {ya.m},  and  all  the  previous  results 
hold  with  respect  to  this  case.  The  question,  however,  is' 
how  to  relate  yn.m  and  y^^  the  latter  being  ultimately 
sought 

This  relationship  depends  on  the  spherical  area  over 
which  the  averaging  is  performed.  It  is  derived  for  the 
uniformly  weighted  angular  block  averages  (e.g.,  l®xl® 
mean  anomalies)  by  (japoschkin  (1980).  Tne  relationship 
for  the  generally  weighted  angular  block  average  is  deduced 
as  follows.  Let  the  average  be  defined  in  terms  of  a 
convolution  with  the  general  weighting  function,  b: 
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b(0,X,e',X‘)  g(0’,X’)  sin©'  d0'  dX* 

.0  Jo 

(51) 

And,  suppose  is  the  frequency  response  of  the 
weighting  function  (the  filter)  in  term  of  the  CarT£5iJ7i 
variables  0,X.  Then  by  the  convolution  theorem  for 
periodic  functions  on  the  plane 

Glc,in  =  Bvijn  Glc,m  (5^) 

Using  (19)  for  both  Gt,m  and  Gk.ni  yields 


g(e,x)=J_ 

4k 


s, 

a=jmi 


:;7 

Yo.ni  Pk 


:jn  ^  Ynjn  Pk 

D=iml 


Qjnil 


(53) 


Now  multiply  both  sides  by  cr!^,  sum  over  k,  and  use  the 
orthogonality  (15)  to  get 


=—  X 

2em 


(54) 


For  purposes  of  discussion,  suppose  the  filter,  b, 
removes  all  power  beyond  the  Nyquist  frequencies  on  the 
plane,  so  that  =  0  for  Ikl  >  K/2  or  Im!  >  M/2.  (The 
response  for  the  uniformly  weighted  block  avrage  is 
rather  inefficient  in  this  respect)  Then,  clearly,  Yj.n  =  0 
for  Im!  >  M/2;  but  ^,m  ^  0  for  any  degree,  j,  when  1ml  < 
M/2.  In  terms  of  the  previous  notation  (see  Figures  1  and 
2),  AYrn  =  0,  but  5Ym  ^  0.  This  means  that  there  remains 
considerable  aliasing  of  the  spectrum  estimated  by  any 
method  of  harmonic  analysis  if  the  averaging  is  performed 
over  constant  angular  blocks.  In  fact,  with  a  perfect 
angular  block  average  (in  addition  to  being  zero  at 
frequencies  higher  than  the  Nyquist,  =  1  for  Ik!  < 
K/2  and  Iml  <  M/2),  equation  (54)  with  (15)  becomes 
Yjjn  ^  Yj.m  for  any  j  if  Iml  <  M/2.  That  is,  the  spherical 
spectrum  denoted  by  the  region  5  y  in  Figure  2  is  passed 
fuDy  by  this  filter  and  therefore  represents  an  aliasing  error 
(see  equation  42). 

A  better  filter  is  the  so-called  spherical  cap  average 
(Jekeli  1981).  In  this  case  (for  constant  size  caps,  only!), 
the  frequency  response  of  the  filter  is  best  formulated  in 
terms  of  its  Legendre  spectrum,  and 

Yn,m  =  Pn  Yn.m  (^^) 


^For  the  uniformly  weighted  average  over  a  spherical  cap  of 
radius, 


b{v)  = 


j  1  -  COSYc 

0; 


otherwise 


(56) 


where  cosvjt  =  cos0  cos0'  +  sin0  sin0’  cos(X  -  X').  The 
filter  response,  is  the  familiar  Pellinen  smoothing 
factor  (Sjdberg  1980).  Figure  3  shows  p^^  for  =  0.5®. 


Jig.3:  Pellinen  smoothing  factor  (thin  line)  versus  Gaussian  smoothing 
factor  (thick  line). 


If  the  filter  were  designed  such  that  p^  =  0  for  n  >JC 
then  from  (55),  clearly  Ya.m  =  0  for  n  >  K  (i.e.,  both  Ay- 
=  0  and  Sym  =  0);  and  the  Legendre  spectrum  of  the 
averaged  function  estimated  by  (31),  (41),  or  (44)  is  not 
aliased.  Such  a  filter  is  appro.ximately  given  by  the 
Gaussian  weighting  function  (Jekeli  1981): 

b(Y)  =  ;  a>0;  0<v<:t  (57) 

The  corresponding  response  is  also  shown  in  Figure  3  for 
a  =  40263.  Both  responses  were  designed  to  attenuate  all 
harmonics  with  n  >  360  to  less  than  20%.  (The  filter 
(51)  with  (57),  being  a  global  average,  is  approximately 
a  cap  average  since  b(Y)=^»  \{;>Yc(a).) 

It  is  seen  from  (54)  and  (55)  that  the  frequency 
responses  of  the  constant  angular  block  average,  B^,  and 
of  the  cap  average,  p^^,  are  completely  different,  the  first 
capable  of  removing  Ay,  only,  and  the  second  able  to 
remove  both  Ay  and  Therefore,  also  the  practice  of 
approximating  the  response  of  the  constant  angular  block 
average  by  pj^  is  quite  inconecL 

Data  Errors  and  Least-Squares  Collocation 

So  far,  the  data  measurement  noise  and  an  a  priori 
spectrum  with  weights  have  not  been  considered.  These 
are  important  aspects  of  modem  harmonic  analyses. 
However,  the  computational  tractability  of  the  analysis 
for  high-degree  models  (such  as  N  =  360)  restricts  the 
allowable  variances  and  covariances  among  the  errors  of 
the  given  data,  as  well  as  the  weights  of  the  a  priori 
harmonic  coefficients  (Colombo  1981).  In  essence,  the 
orthogonality  of  the  matrices  Ym  (see  41)  must  be  availed 
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to  reduce  the  dimensionality  of  the  normal  matrices  to  be 
inverted. 

If  (38)  is  written  as 

g  =  YY  '  (38) 

where  Y  »  Yx-ijandy  =[Y.N+i,  YN’iJ^> 

the  least-squares  solution  with  weights  based  on  the  noise 
variance/covariance  matrix,  P'^,  is  given  by 

Y=(y  P  Yj  Y  Pg  (59) 

In  order  to  take  advantage  of  the  orthogonality  (40),  it  is 
required  that 

(yLpYz)  =  0;  Urn  (60) 

This  happens  only  if  the  errors  in  the  data  {g^  j]  are 
uncorrelated  in  longitude,  and  if  the  correlation  in  latitude 
is  the  same  for  all  longitudes. 

The  absence  of  redundancy  of  data  in  the  new 
approach  obviates  weighting  the  data,  and  no  restrictions 
need  be  imposed  on  the  correlations  in'  the  data  noise  for 
either  formulation,  (44)  or  (31).  The  variance/covariance 
matrix  of  the  estimated  coefficients  is  obtained  by  simple 
error  propagation. 

Colombo’s  (1981)  least-squares  collocaiion  (Ls.c) 
method  of  harmonic  analysis,  in  theory,  is  supposed  to 
yield  the  best  estimator  of  the  Legendre  spectrum  for  a 
given  set  of  data.  In  case  the  data  (again,  globally  and 
regularly  distributed  on  a  grid)  are  perfect  measurements 
(no  measurement  noise),  Ls.c.  minimizes  the  aliasing 
error.  The  Ls.c.  estimator  is  also  of  the  quadratures  t>'pe, 
as  defined  by  Colombo  (1981,  p.43,  eq.2.60): 

=  (61) 
S=0  1=0 

where  i\l^m  are  "weights"  to  be  determined  by  inverting  a 
KxK  matrix  of  discrete  spectral  components  of  the 
Toepliiz/circulant  autocovariance  matrix  of  the  data. 
Although  originally  formulated  to  estimate  only 
coefficients  up  to  degree  and  order  N,  it  is  easily  extended 
to  estimate  any  number  of  coefficients,  such  as  those  in 
(28).  Again,  the  formulation  (61)  is  preserved  in  the 
presence  of  data  noise  only  if  the  noise  covariance  matrix 
has  a  structure  identical  to  the  physical  au  toco  variance 
matrix  of  the  data. 

Considering  (37)  and  (45),  it  is  readily  seen  that  all 
the  other  methods  discussed  here,  namely  (22),  (31),  (41), 
and  (44)  are  of  the  quadratures  type,  but  with  different 
weights;  and  those  in  (61)  are  the  best  in  terms  of 
minimizing  the  aliasing  error  for  a  particular  set  of 


estimated  coefficients.  But  the  optimality  of  (61)  is 
based  on  the  assumed  isotropy  and  homogeneity  of  the 
(physical)  covariance  function  of  the  data.  These 
assumptions  are  quite  strong  and  do  not  hold  uniformlv 
for  all  frequencies.  Furthermore,  the  covariance  function 
usually  is  based  on  a  model  with  questionable  fidelity  at 
the  high  frequencies.  Therefore,  in  practical  applications 
Ls.c.  (for  example,  Colombo  1981;  Gleason  1987)  may 
not  be  significantly  superior  to  the  estimation  based  on 
the  harmonic  analysis  (31)  or  (44)  with  appropriate 
filtering,  e.g.,  with  (57),  because  aliasing  thus  can  be  well 
controlled  and/or  eliminated. 

Summary  ' 

The  methods  of  harmonic  analysis  on  the  sphere  were 
investigated  from  the  point  of  view  of  aliasing.  The 
aliasing  error  was  formulated  in  terms  of  spherical 
harmonic  coefficients  in  the  case  when  the  function  is 
sampled  on  a  regular  grid  of  latitudes  and  longitudes.  This 
yielded  several  results,  some  of  which  are  known,  but 
perhaps  more  analytically  illustrated,  here:  1)  The  simple 
quadratures  method  and  related  methods  are  biased  even 
with  band-limited  functions.  2)  A  new  method  that 
eliminates  this  bias  is  also  superior  to  Colcmbo  s  method 
of  least  squares  in  terms  of  reducing  aliasing,  because  it 
estimates  more  coefficients  at  little  added  computational 
cost.  3)  But,  a  simple  modification  of  the  model  makes 
the  least  squares  method  identical  to  the  new  method  as 
one  is  the  dual  of  the  other.  The  new  method  solves  for 
as  many  coefficients  as  data  and,  in  the  presence  of  data 
noise,  it  need  not  use  fabricated  noise  covariances  to  make 
the  solution  numerically  tractable.  4)  The  essential 
elimination  of  aliasing  can  only  be  effected  ’?.ith  spherical 
cap  averages,  not  with  the  often  used  constant  angular 
block  averages. 

Appendix  -  The  Computation  of 
Analogous  to  (2b) 

=  —  X  PnjnCcosBs)  ,  k  =  0, K-1 


Also  define  the  Discrete  Fourier  Transforms  of  cos0  and 
sin0  on  the  interval  (0,::): 

I,  =  X  y  cos9,  .  It.K  =  he  (A2a) 

K  «o 

le  =  sine,  =  L  (A.2b) 

K  s=o 

The  recursion  formulas  for  the  normalized  associated 
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Legendre  functions  are 

Pn  m(y)  =  Cn-l.m  y  Pn-l,n)(y)  “  dn-2.m  Pj-2.m(y)  .  0  <  m  <  U-2 

(A3a) 

P».n-i(y)  =  V2n+  1  y  PB-t.n-i(y) ,  n  >  1  (A.3b) 

Pn.n(y)  =  V  1  -  r  Ps-l.a-l(y)  .0^2 


Po.o(y)  =  i.  Pi.i(y)  =  V3' 


dD-2jn  - ' 


((2n-l)(2n+l) 

(n  -  m)(n  +  m) 

I  (2n  +  l)(n  +  m  -  IXn  -  m  -1) 
(2n  -  3)(n  -  m)(n  +  m) 


Substituting  (A.3a-d)  into  (A.l)  and  using  the 
convolution  theorem  for  the  product  cose*  Pb.-b(cos05) 
then  yields  recursion  formulas  for  the  Pk'  -  It  =  0,  K-1. 


°n.m  _ 

Pk  “  Cn-ljn 


K-1  ^  ,  ,0 


=  V2n  +  1  £  .  n  ^  1  (A.4b) 

j=0 


=  n>2  (A.4c) 

Y  2n  U 


p°-°  =  5(k,0)  ,  p‘k'=V3^i 
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